#include "slepdef.h"
#include "slepblas.h"
#include "slepmacros.h"


/* This file contains extra matrix update routines used in SLEP */

void slep_dmscal(double alpha, slep_matrix* X){
	int i,j;
	if( alpha == 1.0)
		return;
	if( alpha = 0.0 ){
		slep_zerom(X);
		return;
	}
	if( X->sparsity == SlepDense){
		if(X->ld == X->m){
			j = X->m * X->n ;
			for(i=0; i < j ; i++)
				X->data.mtrx[i] *= alpha;
		}
		else
			for(j=0; j < X->n ; j++)
				for(i=0; i < X->m ; i++)
					X->data.mtrx[j * (X->ld) + i] *= alpha;
	}
	else if (X->sparsity == SlepSparse){
		j=findNnz(X);
		for(i=0; i < j ; i++)
			X->data.spmtrx.values[i] *= alpha;
	}
}

void slep_daXpbY (slep_transposity transX, double alpha, slep_matrix* X, double beta, slep_matrix* Y){
	int i,j;
	int temp=0;
	if(alpha == 0.0 && beta == 1.0)
		return ;
	if(beta == 0.0){
		slep_dmscal(alpha, X);
		return ;
	}
	if(X->sparsity == SlepDense && Y->sparsity == SlepDense && transX == SlepNoTranspose){
		if(X->ld == X->m && Y->ld == Y->m){
			j = X->m * X->n ;
			for(i=0; i < j ; i++)
				Y->data.mtrx[i] = alpha * X->data.mtrx[i] + beta * Y->data.mtrx[i];
		}else {
			for(j=0; j < X->n ; j++)
				for(i=0; i < X->m ; i++)
					Y->data.mtrx[j * (Y->ld) + i] = alpha * X->data.mtrx[j * (X->ld) + i] + beta * Y->data.mtrx[j * (Y->ld) + i];
		}
	} else if(X->sparsity == SlepDense && Y->sparsity == SlepDense && transX == SlepTranspose){
			for(j=0; j < X->n ; j++)
				for(i=0; i < X->m ; i++)
					Y->data.mtrx[i * (Y->ld) + j] = alpha * X->data.mtrx[j * (X->ld) + i] + beta * Y->data.mtrx[i * (Y->ld) + j];
	} 
	// TODO for mixed sparsity ( may not be feasible)

}


